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Variational treatment of the Shastry-Sutherland antiferromagnet using Projected 

Entangled Pair States (PEPS). 

A. Isacsson^ and O. F. Syljuasen^ 

^NORDITA, Blegdamsvej 17, Copenhagen 0, DK-2100, Denmark 
(Dated: February 6, 2008) 

We have applied a variational algorithm based on Projected Entangled Pair States (PEPS) to 
a two dimensional frustrated spin system, the spin-1/2 antiferromagnetic Heisenberg model on the 
Shastry-Sutherland lattice. We use the class of PEPS with internal tensor dimension D — 2, the 
first step beyond product states {D — 1 PEPS). We have found that the D — 2 variational PEPS 
algorithm is able to capture the physics in both the valence-bond crystal and the Neel ordered state. 
Also the spin-textures giving rise to the magnetization plateaus seen in experiments on SrCu2 (603)2 
are well reproduced. This shows that PEPS with the smallest nontrivial internal dimension, D = 2, 
can provide valuable insights into frustrated spin-systems. 

PACS numbers: 



I. INTRODUCTION 

Spins with antiferromagnetic interactions prefer oppo- 
site alignment. However, in many materials the lattice 
structure does not allow all antiferromagnetic bonds to 
be satisfied simultaneously. These are known as frus- 
trated antiferromagnets and display a variety of differ- 
ent phases ranging from rather well-known Neel-ordered 
phases to much less understood exotic phases such as va- 
lence bond crystals and spin liquids. There is no general 
theory of frustrated antiferromagnets, thus the different 
lattice structures are usually studied as separate models. 
Especially difficult are 2D models. Finding the phase 
diagram of any of these models is made difficult by the 
lack of effective numerical tools that goes beyond exact 
diagonalization of the Hamiltonian. 

It is an unfortunate fact that the most powerful numer- 
ical methods such as the Density Matrix Renormalization 
Group (DMRG) and Quantum Monte Carlo (QMC) that 
are very effective for studying general quantum magnets 
do not work well when applied to frustrated 2D antifer- 
romagnets. DMRG is mainly restricted to ID systems, 
and QMC suffers from the sign-problem. However there 
are promising variational mcthodsiiS*^*^ that performs 
an energy minimization in a large class of states known 
as Tensor Product States (TPS). Recently Verstraete and 
Cirac suggested an alternative minimization strategy in 
this space of states, there termed Projected Entangled 
Pair States (PEPS), that promises to be very efficient^ 
and deserves further study. The PEPS or TPS have a 
natural "refinement" parameter, the internal dimension 
D of the tensors. This parameter determines how well 
the particular class of states covers the full Hilbert space. 
The lowest level D = 1 corresponds to product states, 
thus yielding mean field theory results. The aim of the 
present article is to investigate how well the next level 
in the hierarchy, D = 2, can describe a 2D frustrated 
antiferromagnet of real physical interest. 

An interesting frustrated antiferromagnet that exhibits 
both a Neel ordered phase and a valence bond crys- 



tal phase is the spin-1/2 Heisenberg antiferromagnet 
on the Shastry-Sutherland lattice, see Fig. ^ This 
model was initially proposed as a toy model possess- 
ing an exact dimerized eigenstate known as a valence 
bond crystal''. However the interest in this model is 
more than academic as it is believed that the material 
compound SrCu2 (603)2 is reasonably well described by 
this model for particular values of the antiferromagnetic 
couplings^. Although extensively studied, the zero tem- 
perature phase diagram of the Shastry-Sutherland anti- 
ferromagnet remains elusive. While two of the phases are 
known, the possible existence of an intermediate phase 
and its nature are still unresolved issues. In addition ex- 
periments on SrCu2 (603)2 in a magnetic field show the 
appearance of magnetization plateaus^ with rather pe- 
culiar spin structures'^. Several theoretical approaches, 
based on the Shastry-Sutherland model have attempted 
to explain these stepa^i'^i^'^i'^i'^i'^i^'^i'^i". The ap- 
proaches used so far have ranged from exact diagonaliza- 
tioniSiiiiiS, perturbative analysisi2iiiiiiiSiiSiii and mean 
field theory calculations-'^^. 

The PEPS or TPS arc higher dimensional general- 
izations of Matrix Product States^SiSi which are known 
to be particularly useful variational states in 1D2S. In 
contrast to the variational algorithm proposed in Refi^ 
the variat iona l calculations using TPS carried out in 
refs. nUals lilsl build in translational invariance at the 
outset in the minimization procedure by using site- 
independent tensors. While this reduces the number of 
variational parameters it is often desirable not to assume 
this when dealing with a spin system where the a priori 
unknown magnetic unit cell can be bigger than the unit 
cell of the lattice. There are also systems for which trans- 
lational symmetry is explicitly broken by for instance im- 
purities or boundaries. Thus it is desirable to have a 
method that is capable of treating also these situations. 
As a relevant example here, the NMR experiment on the 
1/8 magnetization plateau in SrCu2 (603)2 showed that 
the results were best explained in terms of a state that 
breaks translational symmetrjiiS. 

From the viewpoint of ID variational calculations 
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FIG. 1: Bond configuration for the Shastry-Sutherland model. 
All bonds have antiferromagnetic couplings. Vertical and hor- 
izontal bonds a coupling strength J\ and diagonal bonds J2. 
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FIG. 2: Computational lattice and associated tensors. Shown 
here is an example of a 3 x 3 computational lattice. With each 
site i is associated 2 tensors A^, s =t, j. Each tensor has 
indices u,d,l,r corresponding to bonds connecting the site i 
to neighboring sites. 



where the MPS have internal matrix dimensions D ^ 
32 — 128^^, it would at first sight seem inadequate to re- 
strict the 2D calculations to D = 2. However as argued 
in Ref. |2J even the D = 2 class of states is very rich, a 
fact that is supported by our findings. We find that the 
D = 2 PEPS capture most of the known physics of the 
Shastry-Sutherland model such as the valence bond crys- 
tal phase, the Neel-ordered phase and the magnetization 
steps. 

The outhne of this paper is as follows; In section Hll we 
give a detailed outline of the variational method and in 
Section IIIII we introduce the Shastry-Sutherland model 
and give a brief account of what is known about the 
ground state and the connection to the experimental 
results on SrCu2 (603)2 . In section Hvl we comment 
on the application of variational PEPS to the Shastry- 
Sutherland model and in sections and IVII we look at 
the ground state with and without of external field re- 
spectively, examining the phase transition and magneti- 
zation plateaus. Finally in section IVIII we address the 
performance of the algorithm. 



II. VARIATIONAL METHOD USING PEPS 

Although the algorithm is described in Ref.|6|we reiter- 
ate it here in detail for completeness. As any variational 
algorithm, the aim is to minimize the expectation value of 
the Hamiltonian within a given class of trial states. The 
class of states used here are Projected Entangled Pair 
States (PEPS) represented by an array of complex ten- 
sors Ai, each tensor associated with a physical spin. To 
define a PEPS trial wave function an auxiliary lattice, the 
computational lattice, is introduced. While the sites on 
the computational lattice coincide with the sites on the 
physical lattice the bonds need not. However, it is impor- 
tant that the dimensionality of the computational lattice 



is the same as the dimensionality of the physical lattice. 
The bonds in the computational lattice determine the 
index structure of the tensors Ai. A tensor Ai will have 
one index for each bond in the computational lattice em- 
anating out from site i. Note that different choices of the 
underlying computational lattice lead to different classes 
of variational PEPS-wave functions. As an example, con- 
sider Fig. 12 where a computational lattice in the form of 
a simple 3x3 square lattice with open boundary con- 
ditions is shown. With each lattice site i we then asso- 
ciate two tensors Af, s =t, J. corresponding to spin up 
and spin down respectively. Each tensor has a rank de- 
termined by the number of bonds in the computational 
lattice connecting the site and a dimension D. Hence on 
site 5 in the lattice in Fig. |21 we have two (55 =1"j i) D- 
dimensional rank 4 tensors [^5^];j" with indices for the 
bonds going down, up, left and right, whereas on site 6 
we have two rank 3 tensors [Ag^]f'". While for PEPS one 
associates the tensor indices with the bonds in the com- 
putational lattice one could alternatively associate them 
with plaquettes as in the interaction-round-face TPS^^. 
For a system with M sites we have the following form 
of the trial wave function 



!*)= E ••• E ^^(^1' 
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The symbol Tr{-) means here that one should trace over 
all indices (bonds) in the computational lattice. As an 
example, for the 3x3 lattice in Fig. [5] this operation 
becomes 
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where repeated indices should be summed over. 

For D = 1 the Af:s are complex scalars and the trial 
wave function is a simple product state Ansatz similar to 



a mean field 



c-l c c+1 



M 



!*^=i) -HE A'' i^« 



For D — 2 each index takes on two values. Although we 
will not make explicitly use of it in the following, each 
A^' can for D = 2he represented as a vertex with arrows, 
one for each index, each pointing either in or out. The 
contraction of all indices corresponds then to evaluating 
the partition function of a particular vertex model where 
A*' represent the vertex weights^. 

To minimize the energy (or to even calculate it) we 
need to evaluate 



{H) = 






To see how this is done in practice we consider first the 
normalization A^ = (^|^) with our 3x3 example above 
which explicitly gives 
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We now single out a specific site, say fc = 5, and construct 
the D^ dimensional tensors Ei,i ^^ 5 
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FIG. 3: Steps in obtaining contribution to the effective oper- 
ator matrices A/^ and Ti.% in Eq. (|1J . For a site k located at 
row r and column c we first contract row-wise from the top 
down to row r + 1 and upwards from the bottom to row r — 1 
all the i?-tensors contributing to the operator. After each row 
contraction the left-right dimension of the tensors are reduced 
to a dimension Df before next row is contracted. When only 
rows r and r ± f remain we contract vertically to columns 
c ± 1 after which the effective matrices can be obtained. 



By treating the 2D^ components of ^5 as a 2D^ dimen- 
sional vector A5 and the 4D* components of Ss' ,^5 A^5 as 
a (2D'') X (21)4) inatrix J\ff Eq. ® becomes 

(*|*) = AjAAf A5. 

The evaluation of ($|iJ|^) can be done in a similar 
fashion if we treat each term in the Hamiltonian individ- 
ually, i.e, H = ^^ _ff ("' where 7?^"^ can be written as a 

product of on-site operators i/*^"' — Y[j=i Oj ■ Again 
we form D'^ dimensional tensors i? from the D dimen- 
sional tensors A^^ by 
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Here the tensor product acts on all indices in the tensor, 
i.e., the tensors Ej have composite indices 



For each term _ff *•"■* we can now again write this in vector 
form if we single out a particular site k 
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The normalization can now be written as 
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and sum the matrices "Hj^ to obtain an effective Hamil- 
tonian matrix Ti'{^ — J2n ^i" ^'^^ ^^*^ ^■ 

The overall structure of the optimization algorithm is 
now the following. We pick a site k and calculate J\f^^ 
and Hl!^ by contracting all indices of the i?-tensors sur- 
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Contracting all indices except those connecting site k 
we get 
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rounding it. 
problem 



Then we solve the generalized eigenvalue 



HfA, = AA/-f Afe 



(4) 



from which a new A^ with lower energy can be deter- 
mined. While it is in principle possible to chose this new 
A^, to be the eigenvector corresponding to the smallest 
eigenvalue in Eq. Q) this occasionally leads to problems 
with convergence. Instead we only gradually project out 
the high energy eigenvectors from A^ in the optimiza- 
tion. We then continue in this vein sweeping over all 
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FIG. 4: Boundary conditions used in the simulations, (a) 
Open boundary conditions, (b) Periodic BC. 



sites 1 < k < M until no further reduction in energy can 
be achieved. 

While it is no problem to obtain A/]:* and T-ff.^ in the 
small 3x3 example above it becomes a problem as we 
move to larger systems. For a general set of tensors Ei 
the trace '?"r(J|^ Ei) is in its most general form, an NP- 
complete problem24 and cannot be evaluated exactly for 
large systems but approximate strategies have to be used. 
We have employed the strategy suggested in Ref. 6 doing 
this in a row-wise fashion. To calculate either Af^^ or 

one of the contributions TiJ." the three steps in Fig. Q 
are performed. For a site k located at row r and column 
c we first contract vertically from the top down to row 
r + 1 and upwards from the bottom to row r — 1. When 
only rows r and r ± 1 remain we contract vertically to 
columns c ± 1 after which the effective matrices can be 
obtained. 

When contracting vertically the left-right dimension of 
the i!^-tensors will increase. For instance contracting the 
D^ dimensional tensors [-^2]"^ '^ith [En+2]i '" in Fig. 
generates a new tensor 
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with D"^ dimensional left-right indices. This leads to 
an exponential growth of the left-right index dimensions 
with each row-contraction. To handle this we use the ap- 
proximation technique suggested for approximating an 
MPS with dimension Di with another MPS with a lower 
dimension Df < Di described in Ref. 6, which has also 
been successfully used to simulate time-evolution in ID 
systems^^. Below we give an example of how this is done 
for the contraction of row 1 with row 2 in Fig. |21 

Any row in Fig.|3|can be viewed as a Matrix Product 
Operator (MPO). For the bottom row this corresponds 
formally to a vector 

represented by a set of ND'^ D"^ x D^ matrices E'^. One 
of the middle rows, for instance the second one, can be 
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Contracting row 1 with row 2 is thus formally equiv- 
alent to a vector-matrix multiplication giving rise to a 
new ND^ dimensional vector U21 = U2U1 represented 
by ND^ D" X D" matrices F^ [cnf. Eq. ©]. We now 
seek a new vector 

U2i^Y.^i"---^N"Wir-- ,un) (6) 

represented by ND'^ Df x Df matrices {Df < Z?") F"' 
such that 
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is minimal. One does this in an iterative way starting 
with an Ansatz for the solution and then optimizes the 
matrices f "' one by one until convergence is reached. In 
practice wc do this by first forming the D? x £)? matrices 
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the DfD^ X DfD"^ matrices 

u 

and the D^ x D^ matrices 

u 

in terms of which n can be written 

K = n G. - 2 Re n ^» + n •^- 

i i i 

For a given site 1 < fc < A'^ along the row we can now get 
a linear equation for F1^ that will locally minimize n. To 

see this we differentiate with respect to ( [Ffe]J5 ^, 1 
Ok 
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and treat the left-right indices of F^ and F^ as the in- 
dices of vectors F^ and F^ which leads to the system of 
equations 

GFl = H¥l. 

Although the matrices Ji are not needed to actually do 
the minimization we still calculate them to keep control 



of the error. If the error, after the F^:s have converged 
is too large we increase Df to obtain a better approxi- 
mation. 

For an A^ X A^ system we need to calculate of the or- 
der of N"^ contributions Ti*^"^ to the effective Hamiltonian 
Ti, . For each contribution the contractions and approx- 
imations of MPO:s [cnf. Fig E] need to be calculated. 
This is the most computationally costly part of the algo- 
rithm and to avoid unnecessary calculations we optimize 
the tensors Af in the computational lattice row-wise and 
store all calculated MPO:s which can be reused. This 
means that the memory needed for storage of MPO:s 
scales as N'^D^Dj. 

Finally we would like to point out that in this imple- 
mentation we make no use whatsoever of any symmetries 
of the Hamiltonian, neither in algorithm nor in the trial 
states. This means that our program can treat very gen- 
eral Hamiltonians with nonuniform ground states. 



III. SHASTRY-SUTHERLAND MODEL 

The Shastry-Sutherland model was originally intro- 
duced as an example of a model with an exact dimerized 
ground state'^. The model is a frustrated spin-1/2 anti- 
ferromagnet with a bond configuration shown in Fig. Q] 
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Although the model was introduced for reasons of purely 
theoretical nature, interest was renewed along with ex- 
periments on SrCu2(B03)2^- In SrCu2(B03)2 the crystal 
structure is layered with alternating planes of CuBOa and 
Sr and the magnetic properties stem from the CuBOa- 
layers. It has been argued that these layers are well mod- 
eled by the Shastry-Sutherland mode'^. 

In the limit J2 ^ Ji the ground state, which is sepa- 
rated from the excited states by a gap, is a dimer state 
with localized spin singlets on the diagonal bonds, the 
ground state energy per spin being -Edimor = — 3/8 J2 per 
diagonal bond. In the other limit Ji ^ J2 the model re- 
verts to the ordinary antiferromagnetic Heisenberg model 
with a Neel ordered ground state and gapless spectrum. 
From high temperature series expansion and exact diago- 
nalizationpiii2Si2£ a possible direct transition between the 
dimer phase to the Neel phase has been estimated to lie 
at (Ji/J2)c = 0.7 ±0.01. 

Other works point to the existence of an intermediate 
phase between the antiferromagnet and the dimer phase. 
A sketch of the phase diagram is shown in Fig. El Most 
estimates agree that below (Ji/J2)ci > 0.6 the ground 
state is the dimer state and above (Ji/J2)c2 < 0.9 the 
ground state is the Neel state. The nature of this inter- 
mediate state has been addressed in several publications. 
In Ref . W Albrecht and Mila used Schwinger Boson mean 
field theory to argue in favor of a first order transition 
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FIG. 5: Proposed phase diagrams for the Shastry-Sutherland 
model. Two possible scenarios for the phase diagram have 
been proposed. Either a direct transition between the dimer 
state and the Neel state (a) or a transition via an intermediate 
phase as in (b). The nature of this intermediate phase has not 
been established. 



between the dimer state into a helical state and a sec- 
ond order transition to a Neel state. Another possible 
intermediate state, the plaquette singlet phase was dis- 
cussed by Koga and Kawakimi in Ref. 31. Both plaquette 
states and Helical states were considered in Ref. [3^ Ar- 
guments against bot h p laquette and helical phases was 
put forward in Ref. 33 who performed extensive series 
expansions around both the helical and plaquette phases 
and even columnar phases. This is in contrast to Ref. l34l 
which supports either a plaquette phase or a columnar 
phase. Other suggestions for the intermediate phase are 
Weakly Incommensurate Spin-Density Waves (WISDW) 
or Fractionalized Quantum Para-Magnet (FQPM)^. Fi- 
nally a resonant valence bond plaquette phase was sug- 
gested as the intermediate state in Ref. "Sff. Thus, neither 
the existence of an intermediate phase nor its exact na- 
ture are presently known. 

Experiments on SrCu2(B03)2 in strong external fields 
show the existence of magnetization plateaus^. While 
the ground state in absence of an external field is believed 
to be the dimcr-state, the magnetization steps were orig- 
inally thought to be formed by strongly localized triplets 
forming a periodic patterns which may spontaneously 



break the translational symmetrj* 



12.13.14.15.16.17 



(for an 



alternative explanation hypothesis see Ref. IT8|) . Subse- 
quent NMR experiments^- at the 1/8 plateau revealed 
a more complex structure, inconsistent with the simple 
triplet-singlet picture. By including coupling to phonon 
degrees, with the sole purpose of breaking the transla- 
tional symmetry, exact diagonalization studies of small 
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revealed more complex spin textures. Again, 



at all steps (except the 1/2) translational symmetry is 
broken and larger unit cells are formed. 



IV. APPLICATION OF VARIATIONAL PEPS 
TO THE SHASTRY-SUTHERLAND MODEL 

Applying variational PEPS to the Shastry-Sutherland 
model is straight forward. Two issues should be noted. 
First, we stress again that the computational lattice does 
not need to have the same bond-configuration as the un- 
derlying Hamiltonian. For the purpose of studying the 
Shastry-Sutherland model it is sufficient to use an ordi- 
nary square lattice as depicted in Fig. [3 It is easy to 
show that already with a low tensor dimension D = 2 
it is possible to represent exactly the dimerized ground 
state with singlets on all diagonal bonds. 

The second issue regards boundary conditions. In 
implementing the algorithm we have used a computa- 
tional lattice with open boundary conditions. The rea- 
son for this is two-fold. Firstly, using a computational 
lattice with periodic boundary conditions severely re- 
duces the performance of the algorithm, the difference 
being that between matrix-vector multiplications rather 
than matrix-matrix multiplications. Secondly, our pro- 
gram suffers from stability problems arising due to ill- 
conditioning and round-off errors in the case of compu- 
tational lattices with periodic boundary conditions. 

Although, the computational lattice docs not have pe- 
riodic boundary conditions it is not necessary to adopt 
the same boundary conditions to the Hamiltonian. In 
this study we have used two different physical bound- 
ary conditions, open and periodic. For the open BC we 
have adopted the geometry shown in Fig. 0fa) while for 
the periodic the geometry in Fig. E^b) . The interpreta- 
tion of using periodic BC in the physical problem while 
using open BC in the computational lattice is reminis- 
cent of using a self-consistent field on the boundary. The 
bonds across the boundary have only a tensor dimension 
^boundary = 1 which in thc limit of large lattices implies 
that the contribution from a such a bond will approach 
the product form (Si • Sj) — > (Si) ■ (S) for limited D. 

We have mainly restricted ourselves to using a tensor- 
dimension D = 2 for which a usual work-station with 
1GB of internal memory suffices. Although our program 
can in principle handle D > 2 (see section. IVIip . it is in 
its present incarnation too slow and unstable for D > 2. 
For the dimension of effective MPO:s when calculating 
effective operators we have used a variable 16 < Df < 24 
for all simulations except for the largest (12 x 12) systems 
where the 1GB memory limit restricts us to 16 < Df < 
18. 



V. GROUND STATE IN ZERO FIELD. 

In Fig.inithe lowest energies obtained by the algorithm 
for D = 2 are shown for system sizes of 6 x 6, 8x8, 10x10 
and 12 x 12 with both open (main panel) and periodic 
boundary conditions (inset). The total energy has been 
scaled by the number of internal diagonal bonds N^ and 
the coupling energy J2- As can be seen, for low J1/J2 




FIG. 6: (Color online) Variational minimum energy as a func- 
tion of J1/J2 for system sizes 6x6 (blue squares), 8x8 (red 
circles), 10 x 10 (green diamonds) and 12 x 12 (black stars). 
For open boundary conditions (main figure) the transition 
from the dimer state to the Neel state is clearly visible in 
the energy which has been scaled to the number of diago- 
nal bonds. The inset shows energies for periodic boundary 
conditions using the same energy scaling. 



the ground state energy approaches -3/8J2 for the sys- 
tem with open boundary conditions confirming the con- 
vergence to the dimerized ground state. From the graph 
obtained using open boundary conditions it can be seen 
that already for 13 = 2 we find a phase transition from 
the dimerized state. The transition point being located 
at 0.69 ± 0.02 which is in agreement with estimates for 
the transition point of the direct dimer-Neel transition. 
Further, in the energy a clear finite size effect is seen as 
the transition point is approached. 

From the inset showing the energies obtained using 
periodic boundary conditions the location and nature 
of the transition is less clear. Here we have again di- 
vided the total energy by the number of internal diago- 
nal bonds. Since the diagonal bonds across the bound- 
ary are not counted this gives rise to energies lower than 
— 3/8J2. The transition is revealed by looking at the 
strength of the diagonal singlets and the staggered mag- 
netization, as shown in Fig. [7| The singlet mixing is 
calculated by projecting the diagonal bonds on to the 
singlet state, i.e. 1 indicates a singlet while indicates 
a triplet. The staggered magnetization displayed is cal- 
culated as 2< {Y.,Si{-iyf >i/2. Note that for large 
J1/J2 the staggered magnetization is higher than the 
value expected for a Heisenberg antiferromagnet. The 
reason for this can be two-fold. Firstly, finite size effects 
explain a part of the discrepancy as can be seen from the 
graph. Secondly, while D = 2 gives a good value for the 
energies involved, being only a few percent off the exact 
values, observables may differ by more (see Sec. lVII|) . 

As stated in Section UTTl there are good reasons to be- 
lieve that an intermediate phase exists between the dimer 
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FIG. 7: (Color online) Order parameters, singlet mixing 
on diagonal bonds (solid lines) and staggered magnetization 
(dashed lines) for periodic boundary conditions for system 
sizes 6x6 (blue squares), 8x8 (red circles), 10x10 (green dia- 
monds) and 12x12 (black stars). Results obtained with in- 
ternal tensor dimension D — 2 suggests a direct first order 
transition from the dimer-state to the Neel state. 



phase and the Neel phase. However, it is hardly sur- 
prising that we don't see this intermediate phase in our 
D = 2 variational calculation. First of all, with ten- 
sor dimension D = 2 we typically overshoot the true 
ground state energy by a few percent, thus higher D is 
likely needed to capture any additional phase that may 
differ by a percent or less in energy. Second, the influ- 
ence of boundary conditions scales as 1/A'' which implies 
that boundary effects can have a big impact even for the 
largest systems (12 x 12) in cases where the energy split- 
tings between ground state candidates are small. 



VI. MAGNETIZATION PLATEAUS 

To study the magnetization curve we have deliberately 
chosen a somewhat smaller coupling constant J1/J2 = 
0.6 than the experimental value [Ji/ J2)srCu2{B03)2 — 
0.635. This makes the dimer state more stable and the 
algorithm converges faster but should not significantly 
affect the physics. In Fig. |S| we show the results for a 
finite magnetic field for systems with periodic boundary 
conditions of sizes 8x8 and 10 x 10 (£> = 2). Although 
our square geometry and (periodic) boundary conditions 
are inconsistent with the unit cells proposed in Refs. llOl 
1!? a clear step like structure is nevertheless visible in 
Fig. |S1 A closer inspection reveals that locally the spin 
configurations we have obtained match those in Ref. Il9| 
very well. 

To study the spin configurations at the plateaus we 
have visualized the wave functions by coloring the bonds 
according to the amount of triplet or singlet mixing (see 
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FIG. 8; (Color online) Magnetization curve at J\/ J2 = 0-6 for 
system sizes 8x8 (blue squares) and 10 x 10 (red circles). A 
clear step like structure corresponding to fillings 1/4,1/3 and 
1/2 are seen. The insets show the distribution of singlets (red 
lines), triplets (blue lines). The Sz component of the spin 
on each site illustrated by a circle with radius proportional 
to I (Sz) |. Filled circles are aligned with the field while open 
circles represent spins pointing opposite to the field. 



Fig. ISJ. Localized singlets are drawn in red whereas 
triplets are blue. To determine the color of a bond we 
use the following criterion 




I*) 



< O.letriplet 
> 0.9etriplct 

otherwise 



0.9e; 
O.les 



singlet 
elet 



(8) 



Furthermore we have measured the spin component par- 
allel to the direction of the applied field (B — Bz) and vi- 
sualized {Sz) by circles with radii proportional to | (Sz) \- 
Spins aligned (anti-aligned) with the field are drawn as 
filled (open) circles. 

For small fields, i? < 1, we see a finite magnetization 
where one would expect a spin gap. This is due to the 
inability of our trial states to form singlets across the 
boundary. As can be seen in the inset, for B < 1 the 
interior of the system is still in the dimer phase while only 
spins on the boundary have aligned with the field. Thus, 
by looking at when the magnetization in the interior of 
the system becomes finite we estimate the spin gap to be 
roughly B = 1 corresponding to 33 T where we have used 
coupling constants J = 85 K, g = 2.28 (with J=71K the 
corresponding number is 28 T). 

For B > 1 three steps can be distinguished, 1/4, 1/3 
and 1/2. While the spin texture at 1/2 matches that 
of earlier predictions, half of the diagonal bonds being 
triplets while the other half being singlets, the spin tex- 
tures for other points are more elaborate and only agrees 
with earlier predictions locally. An example is shown in 
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FIG. 9; Convergence of energy in a 10 x 10 system with peri- 
odic boundary conditions for D = 2. The energy is shown as 
a function of the number of diagonalizations of the general- 
ized eigenvalue problem in Eq. |1] The inset shows a closeup 
of the final part of the optimization where fluctuations due 
to the approximation strategy used to calculate expectation 
values are visible. 



the top left inset of Fig. |H1 Here, three of the 1/3 unit 
cells obtained by Miyahara et aU^ are reproduced in the 
interior of the 10 x 10 system. Note that the total mag- 
netization is larger than 1/3 at this point due to the spin 
configurations on the boundary. 

The spin textures obtained at the steps can only be 
found when translational symmetry is broken. In the 
variational method employed here translational symme- 
try is broken partly because of our choice of computa- 
tional lattice and boundary conditions, and partly be- 
cause a D = 2 PEPS cannot represent the coherent 
superposition of degenerate plateau states connected by 
global symmetry transformations. 



VII. ALGORITHM PERFORMANCE 

So far we have only concerned ourselves with D — 2 
which is the first step beyond a simple on-site factorizable 
wave function. Restricting ourselves to D ^ 2 and sizes 
up to 12 X 12 allows the program to run on an ordinary 
workstation with 1GB internal memory without using 
any swapping to disk. 

Because the method is based on a sequence of approx- 
imations, i.e., for a 12x12 system there are over 700 con- 
tributions to the effective Hamiltonian on any given site, 
each contribution being obtained in a series of up to 10 
consecutive approximations one has to ask whether or not 
the precision is compromised. Another important factor 
to consider is how much more accuracy (how much closer 
to the true ground state energy we can come) can be ob- 
tained by increasing D, and how the computational effort 
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FIG. 10: Deviation of the norm from unity during the opti- 
mization of a 10 X 10 system with Ji = J2. 



scales with increasing D. 

In Fig. Othe energy as a function of number of opti- 
mizations is shown for a 10 x 10 system with Ji = J2 = 1. 
While smooth on a large scale, the errors accumulated in 
the successive approximations are clearly visible in the 
inset which shows a closeup of the final convergence. 
In this simulation and others we have used a final di- 
mension in the approximation of S-tensors (see Sec. ^ 
16 < Df < 24 (For 12 x 12 we have been restricted 
to 16 < Df < 18 due to the limited memory (1GB) of 
the workstation). As can be seen, despite the heavy re- 
duction of the state space in the calculations, we have a 
precision of the order of 4 digits. This can also be seen 
by looking at the norm of the wave function. In Fig. 1101 
the deviation of the norm from the nominal value 1 is 
shown for the simulation in Fig. |5| 

To estimate the accuracy, i.e. how close to the true 
ground state energy the variational PEPS algorithm can 
get we have compared it to QMC on an ordinary Heisen- 
berg antiferromagnet (QMC being unable to handle the 
SS- model), H ~ 1"^,^ -K Si- Sj with open boundary con- 
ditions. The QMC was run at an inverse temperature of 
/3 = 256/ J and the results are shown in Fig.llllwhere the 
quantity 1 — iJpEPs/^'QMC is shown for different system 
sizes and tensor dimensions D. For the smallest system 
4x4 sites we find good agreement with the figures re- 
ported in Ref.|y obtained using imaginary time evolution. 
We further note that as the system size is increased the 
relative error decreases slightly, and that a linear increase 
in D seems to give an exponential increase in accuracy. 

We want to point out that although the energies ob- 
tained by PEPS are in good agreement with exact results, 
observables may deviate more. In the left inset of Fig.lTTl 
the staggered magnetization for the Heisenberg model is 
compared with the staggered magnetization obtained us- 



-D=2 
-D=3 





l-e-QMC 1 
\ — D=2PEPS 

V 


4 6 8 10 12 



9 10 11 12 



FIG. 11: (Color online) Comparison between ground state 
energies obtained by variational PEPS and Quantum Monte 
Carlo for open Heisenberg antiferromagnet with open bound- 
ary conditions for system sizes 4 x 4 to 12 x 12 and ten- 
sor dimensions ranging from D = 1 to _D = 3. Inset (a) 
shows a comparison between QMC and D = 2 PEPS for 
the staggered magnetization for the Heisenberg model. Inset 
(b) shows a comparison between exact diagonalization of the 
Shastry-Sutherland model with periodic boundary conditions 
and PEPS _D = 1 to D = 3. 



ing D — 2 PEPS. Compared to the accuracy in energy 
which is of the order one percent (see main panel) the 
error in the staggered magnetization is an order of mag- 
nitude larger. 

Although the above comparison for the Heisenberg 
model does not, in a strict sense, tell us anything about 
the accuracy obtained for the SS-model away from the 
limit J2/J1 » 1 it still serves as a good indication on 
the general behavior as one varies N and D. For the 
SS-model we have compared with exact diagonalization 
results obtained using SPINPACK2LS at Ji = J2 using 
periodic boundary conditions. The comparison is shown 
in the right inset of Fig. ITTI 

The internal tensor dimension D plays an important 
role in how faithfully a PEPS can represent ground states. 
The algorithm scales very badly with increasing D, the 
bottle neck being the contraction of two rows (See Eq.EJ. 
To form an MPO by contracting two rows requires of the 



order ND^Dj {Dj > D"^) and scaling with D is at best 
D^^ . This then sets a limit to the maximum internal 
dimensions that can be practically used and the large 
values of -D ~ 10^ used in ID variational MPS cannot be 
reached. However, as we have seen here (see also R,ef . V2M . 
for 2D PEPS we expect small D ^ 2—5 to be able capture 
the essential physics for many problems with short range 
interactions. The scaling of the algorithm with linear 
system size is N^ and is less severe. 

All simulations were run on Linux workstations with 
2.0 GHz AMD Athlon processor and 1GB of internal 
memory. On such a machine the graph in Fig. El took 
96h to produce. 

Finally we comment on the stability of the algorithm. 
We have found that the technique used to optimize ten- 
sors one by one often becomes unstable and may not be 
the optimal way to find the ground state. It may well be 
that using imaginary time evolution is a more effective 
way. 



VIII. CONCLUSIONS 

We have applied a variational procedure based on Pro- 
jected Entangled Pair States (PEPS) to study the ground 
state properties of a frustrated spin system, the Shastry- 
Sutherland model. Using the smallest nontrivial dimen- 
sion on the tensors I? = 2 a direct phase transition be- 
tween the dimer state and the Neel state can be observed, 
the location being well in agreement with other theoreti- 
cal estimates for a direct transition. Within D — 2we see 
no clear indication of an intermediate phase which may 
require higher D, larger systems, or proper handling of 
periodic boundary conditions. We also find that already 
with PEPS D = 2, magnetization plateaus are possible 
to reproduce, and that the non-trivial spin textures as- 
sociated with these plateaus can be seen. 

Furthermore we have examined the performance of the 
algorithm and conclude that it degrades rather severely 
for intermediate to large values of the internal dimension 
D. However, this scaling of the performance degradation 
might not be so restrictive as already the D = 2 class of 
PEPS is well suited for studies of frustrated spin systems 
at a level beyond mean field theory. The method can 
readily be extended to other 2D frustrated spin-models. 
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